###################
# Post Estimation

library(rstan)
library(foreign)


# Load Datasets and combine in Data.frame
#########

  chc <- readRDS("CHES_POLICY_COMBINED_DATA_06_14_min3perc.Rds")
  nrow(unique(chc[,c("pid","year")]))
  
  c <- unique(chc$ctry_name)
  y <- c("2006","2010","2014")
  
  cd <- data.frame()
  
  for(c in unique(chc$ctry_name)){
    for(y in c("2006","2010","2014")){
      fit <- readRDS(paste("res_",c,"_",y,".Rds",sep=""))
      d <- summary(fit,probs = c(0.01,0.025, 0.25, 0.50, 0.75, 0.975,0.99))$summary
      d <- as.data.frame(d)
      d$par <- rownames(d)
      rownames(d) <- NULL
      d$ctry_nam <- c
      d$year <- y
      cd <- rbind(cd,d)
  }}

# Check for Convergence
######
  
  # Check for Rhat
  sel_case <- unique(cd[which(cd$Rhat>=1.1),c("ctry_nam","year")]) 
  print(sel_case)  # no cases above 1.1 Rhat
  summary(cd$Rhat)
  
  # Check for n_eff
  eff_cases_below <- unique(cd[which(cd$n_eff<400),c("ctry_nam","year")])
  print(eff_cases_below) # no cases below 400 effective draws
  summary(cd$n_eff)
  
  
# Create Dataset with party estimates
#######

  s <- grep("x\\[[0-9]+\\]|nu\\[[0-9]+\\]",cd$par)
  res_x_nu <- cd[s,]
  
  # Put PID
  res_x_nu$pid <- NA
  for(c in unique(chc$ctry_name)){
    for(y in c("2006","2010","2014")){
        dsel <- apply(chc[,paste("pos.",1:17,sep="")],1,function(x) !all(is.na(x)))
        upids <- unique(chc$pid[dsel & chc$ctry_nam==c & chc$year==y])
        res_x_nu$pid[res_x_nu$year==y & res_x_nu$ctry_nam == c] <- rep(upids,2)
        if(length(rep(upids,2)) != length(res_x_nu$pid[res_x_nu$year==y & res_x_nu$ctry_nam == c])) cat(c,y,"\n")
    }}
  
  # Merge party names
  chall <- read.csv("1999-2014_CHES_dataset_means.csv")
  pid_nam_dta <- chall[,c("party_id","party","year")]; names(pid_nam_dta) <- c("pid","pid_nam","year")
  pid_nam_dta$pid_nam[pid_nam_dta$pid==2406] <- "TB-LNNK"
  pid_nam_dta$pid_nam <- as.character(pid_nam_dta$pid_nam)
  
  res_x_nu <- merge(res_x_nu, pid_nam_dta)


  
  # Save  
  saveRDS(res_x_nu,"res_party.RDS")

  
# Create dataset with item estimates
####### 
  
  s1 <- grep("[a-b]{1}k\\[[0-9]+\\]",cd$par)

  res_items <- cd[c(s1),]
  res_items$par
  
  saveRDS(res_items,"res_items.RDS")
 

